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I. INTRODUCTION 

Experimental investigation of solids is in most 
cases concerned with observation of static or dy- 
namic properties in a weakly perturbed macro- 
scopic system. Therefore, standard techniques 
from equilibrium statistical mechanics are usu- 
ally sufficient, possibly supplemented by linear- 
response theories to account for transport. Equi- 
librium statistical mechanics is based on the Gibb- 
sian approach where the statistical density ma- 
trix of a state at energy E s is given by the Boltz- 
mann factor e~^^ H ~^ N ^ with inverse temperature 
j3 = l/k-oT and the chemical potential fi. The big 
success in the theoretical description of quantum 
systems in thermal equilibrium is based on the 
fact that both the thermal average and time evo- 
lution are based on the same operator, and one 
can use the concept of Wick rotation to formu- 
late a theory which actually condenses both type 
of dynamics into a single complex Matsubara fre- 
quency theory. 

The advances in experimental methods over the 
past two decades have however opened the access 
to studies, where time dependencies on the scale 
of internal time-scales become visible,— or where 
mesoscopic systems can be driven out of ther- 
mal equilibrium in a controlled way and various 
properties can be experimentally observed, both 
in steady- and time-dependent states. Therefore, 
one pressing question to modern quantum many- 
body theory is how one can describe generic non- 
equilibrium situations in macroscopic or meso- 
scopic systems. For the latter the paradigms are 
the single-electron quantum dot and nano-wires, 
where a tremendous amount of data on transport 
or transient response has been collected over the 



past ten years 

Out-of-cquilibrium many-body theory is an 
emerging field which poses an extreme challenge. 
There are many attempts to use existing theo- 
retical approaches, the most popular being the 
ones based on the Keldysh formulation of pertur- 
bation theory^ In particular, the growing inter- 
est in transport through mesoscopic systems trig- 
gered a variety of applications of this technique; 
for example direct perturbation theory with re- 
spect to different zeroth order Hamiltoniansr r — 
functional rcnormalization group methods or their 
derivatives^— or direct numerical evaluation of 
the real-time propagatorsAi~— There are many 
other ideas, for example based on the concept 
of infinitesimal unitary transformations^ A com- 
prehensive overview can for example be found in 
Ref. [l7|. 

An early attempt to formulate an out-of- 
equilibrium version of statistical mechanics for 
steady-state properties of general quantum many- 
body systems is due to Zubarevp£ who tried to 
construct a time-independent density matrix for- 
malism by solving the equation of motion within 
the scattering state formalism. This approach 
has later been revisited by Hershficld in the con- 
text of transport through quantum dot systems^ 
The main problem with these, in principle ex- 
act formulations, is that they cannot be readily 
applied because they require the solution of the 
Lippmann-Schwingcr equation for the scattering 
states, which amounts to knowing the full solu- 
tion itself. There have been attempts to tackle 
this problem by utilizing advanced tools of quan- 
tum many-body theory like Bethe ansata^ or an 
extension of Wilson's numerical rcnormalization 
techniques »2i However, the former approach could 
only be applied to a very specific model, while the 
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latter may lack a thorough foundation regarding 
the proper steady-state limit ^ 

In the present manuscript we focus on a dif- 
ferent way to extend the theoretical framework 
of equilibrium quantum mechanics to steady-state 
nonequilibrium for quantum impurity models via 
an imaginary-time theory. We especially dis- 
cuss the possibility to deform the complex time 
contour for physical observables in equilibrium 
to the Keldysh contour appropriate for nonequi- 
librium, as proposed by Doyon and Andrei^. 
One fundamental problem that arises in any 
such attempt stems from the fact that the non- 
equilibrium steady-state Boltzmann factor and 
the time-evolution operator now have a funda- 
mentally different structure, and thus a straight- 
forward Wick rotation is not possible. As an alter- 
native procedure, we show that, by introduction 
of Matsubara voltage, the problem of the dual op- 
erators can be resolved and a consistent theory 
for steady-state non equilibrium based on auxil- 
iary statistical mechanical problems formulated. 

As the first step we need to properly define in 
what sense we achieve a steady state in a quan- 
tum impurity model. This is done in section [TT] 
together with a discussion of the general struc- 
ture for Keldysh perturbation theory, the prob- 
lem of analytical continuation and the idea of the 
Matsubara voltage formulation. The equivalence 
of the Keldysh real-time and the Matsubata volt- 
age perturbation theory for the steady state will 
be shown in section [TTTl for the single-impurity An- 
derson model. In section ITVl we derive expressions 
for calculating static observables on the impurity 
via an analytical continuation procedure from the 
Matsubara voltage description. As summary, sec- 
tion |V] will conclude the paper. 

Since many details are rather technical and not 
really necessary to understand the main line of ar- 
gument, we included them in a series appendices, 
which will be referred to when necessary. 



II. MANY-BODY THEORY OFF 
EQUILIBRIUM 

A. Convergence to steady-state 
nonequilibrium 

To establish a steady-state nonequilibrium, one 
requires the system to be in the infinite-size limit. 
In mesoscopic systems, such as quantum dots, 
this requirement means that the size of the reser- 
voirs, L, should be the largest scale. In partic- 
ular, the time tw for the wake of the perturba- 
tion occurring in the quantum dot region to reach 



the edge of the reservoir with the Fermi veloc- 
ity Vf (tw = L/vp) should be greater than any 
time scale used for the turn-on of the perturba- 
tion or measurements. This ensures that the re- 
flected wave does not interfere with the formation 
of the steady-state and its measurements. Alter- 
natively, the reciprocal vf/L also represents the 
level spacing of the continuum states, which sets 
the smallest energy scale in the model. 

As in conventional many-body theory, we con- 
sider a perturbation which we turn on infinitcsi- 
mally slow with a rate r/ _1 as 

V(t) = Ve vt (1) 

for the time interval t G [— T, 0], where T is some 
initial time which eventually will be sent to infin- 
ity. For t > 0, the perturbation remains constant 
at the full strength, V(t) = V. The above discus- 
sions lead to the relation between the three energy 
scales 



In his original proposal^, Hershfield assumed the 
presence of an external relaxation process to de- 
rive the time-independent density matrix in the 
limit T —> oo. Recently Doyon and Andrei^ 
have shown that for mesoscopic systems infinite 
reservoirs provide a relaxation process and any 
assumption of an additional external relaxation 
source is not necessary. This suggests that we can 
do away with the adiabatic factor e* 7 * in a time- 
dependent theory as long as the limit L — > oo is 
taken first. Here we show through an explicit cal- 
culation that the adiabatic factor e nt is not nec- 
essary for the steady-state if local measurements 
are made near the quantum dot^l, henceforth ab- 
breviated as QD. 

Our model system consists of a QD connected 
to two fermionic reservoirs labeled by a — L, R (or 
±1, respectively, when the reservoir index is taken 
numerically). We include the single- particle tun- 
neling between the leads and the QD into the non- 
interacting part of the Hamiltonian, which then 
becomes the resonant level model (RLM) 

H = ^ ^akaC^Caka + (d ^ d ° d ° 
aha <J 

~^2~^^ Caka + h - c ■ ( 3 ) 

aka 

Here, c ak is the creation operator of conduction 
electrons for the reservoir a with energy e a ka at 
the continuum index k and spin a; creates an 
electron on the QD orbital and t a is the tunneling 
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integral. Q is the normalization due to the volume 
of the reservoirs. This Hamiltonian can be diag- 
onalized by the scattering state operators ^ aka 
given by the formal Lippmann-Schwinger opera- 
tor equation, 



1 



4, (4) 



with the Liouville operator acting on the operator 
space as CqO = [Hq,0]. This equation can be 
easily solved as 



rex. 



t 



^aka 



4f n e aka - - e a >k'cr + z0+ 
a' h' a 



' (5) 



with the bare retarded Green's function for the 
QD, go(u) = (w - t d + tT)- 1 . Here, V = n(t 2 L + 
t 2 R )N(0) is the hybridization broadening, and we 
assume for simplicity a flat density of states N (0) 
for both reservoirs. 

According to Hershfield^, the nonequilibrium 
steady-state created by a shift of chemical poten- 
tial on the source (drain) reservoir by $/2 (— $/2) 
can be described by a density matrix 



Po 



exp[-P(H - $%)] 



with the so-called F-opcrator defined as 



(6) 



(7) 



and the generalized partition function 
S = Trcxph/3( J ff -*r )] . 

Since Yo is diagonal in the eigen-operators, 
[Ho,Yo] = and po is time- independent. It is im- 
portant to realize that the convergence factor i0 + 
in the denominator of the Lippmann-Schwinger 
equation determines that the one particle states 
c f 

u afcer 



originate from the infinite past inside the 
reservoir of infinite size. Thus the limit L — > oo 
has already been taken implicitly before the per- 
turbation is turned on. 



B. Real-time theory for open system 

In addition to the nonintcracting part Hq, the 
full Hamiltonian H of the system will in general 
also contain an interaction we will denote as V 



in the following. For a general observable A, we 
define its nonequilibrium expectation value as 

Tr (e i6T Ae- i6T po) 
lim i — J -, (8) 



lim (A(T)) 

T— >oo 



T->oo 



Trpo 



where A has been evolved with the full Hamilto- 
nian H during the time interval — T < t < 0. Un- 
like Eq. Q, here we take V{t) = V for —T <t< 
0. Defining the time-dependent operator A(t) in 
the Heisenberg picture, A(t) = e i6t Ae~ i6t , A(t) 
satisfies j- t A(t) = i[H,A(t)] and 



A(t) = A + i / dt'[H,A(t')]. 



(9) 



One can now form the average with respect to po, 
to obtain 

(i(T)) = (A) +i [ dt'([H,A(t')]} 

J-T 

= (A) + i f dt'([v,A(t')}) . (io) 

J-T 

For the existence of a well-defined limit (A(oo)), 
one must show that^ 

fO 

dt([V,A(t)}) <+oo . (11) 

To this end one argues that as long as V and A are 
operators local to the quantum dot,— the time- 
evolution of A(t) will decay as electrons travel 
away and the integral is finite. 

To make the argument concrete, we consider 
as example the usual on-site Coulomb interaction 
V = Und-fUdi and measure the current through 
the dot, A = I. The occupation number operator 
can be expressed in terms of i^aka as 



kk'aot' 



-9d( e k)9d{e'k)^l k<J ^a'k'a- (12) 



With the requirement that the current through 
the L/R leads, Il/Ri is the same, the current 
operator I can be symmetrized as I = (t 2 R lL + 
arid 

_ —ithtR 



tht R i 

L R kk 



tLCRka)) - h.C.] (13) 

Y,(9m-uh')) 



tLtR{l/>l k 1p 



>Lk' 



tpRk^Rk' 



"(*£ - t 2 R )(lpi k 4>Rk> + iik^Lk')} (14) 
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We evaluate Eq. (fTU)) using Wick's theorem. Due 
to the commutator inside the expectation value, 
only connected contractions between any V and 
I(t) will contribute. Therefore any non-vanishing 
Wick's contractions must have an even number of 
contractions connecting V and I(t) and contain 
a factor (ip a ka(0)^l klT {t))o or (ipl krT {0H' aka (t)) o . 
More specifically, the first order perturbation in- 
volves factors like 

([V,Io(t)]) ex ^5>2(fc) - 9 d (k'))g d (k)g* d (k' 



■■(fLk-fLk>)e- l ^-^ t + --- 



(15) 



Summation over the continuum variables k, k! 
leads to terms of the form 



1 



(it(f)d(O)) = -Y^9 d {k)Uk)e-^ 

k 



-ie k t 



k 



(16) 



Note that the inequality holds both for equilib- 
rium and nonequilibrium. Therefore, the follow- 
ing expression 



([V(s k ),[--- ,[y( Sl ),j (t)]---]>o 

„ -r-min{\ Sl -t\,- ,\s k -t\} 



(17) 



holds to any order of the pcrturbative expansion 
in V, and the integral over t, Eq. (fTTj). becomes 
convergent. This shows that the steady-state limit 
of the nonequilibrium is well-defined regardless of 
the adiabatic rate r\. 

However, it should be emphasized that, al- 
though the convergence factor e vt is not necessary 
for a time-dependent theory, such adiabatic factor 
should be treated carefully in a time-independent 
theory, like the steady-state nonequilibrium. Such 
situation arises in particular when we perform 
a Fourier transformation to represent a steady- 
state quantity in a spectral representation with 
sinusoidal basis. For instance, let us express a 
steady-state quantity A as an integral over a time- 
dependent function F(t), 



A 



F(t)dt, 



(18) 



where the integral is absolutely convergent with- 
out any adiabatic factor e nt . We write F(t) in a 
spectral representation as 



(19) 



(a) 



t = -T 



t = -T-ip 



(b) 



t = 

t = -i/3 



t = 

t = -i/3 



FIG. 1: (a) Keldysh contour for real-time di- 
agrammatics. If the time-evolution along the 
dashed line does not contribute an extra factor, 
the whole contour can be deformed to one along 
the imaginary-time from t = to t = as 
shown in (b). 



with the Fourier component F(u>), and the quan- 
tity A becomes 



r° , r r duj ~, , 

A= dt -F(w)e- 

J— oo J— oo ^ 



(20) 



If we now want to express A via a spectral rep- 
resentation, we need to change the order of in- 
tegrals. However, e~ lujt is an oscillatory function 
and we have to insert a regularization factor e vt to 
unambiguously allow the integral exchange. Then 



.4 



did ~, . 
Z7T 

du iF(ui) 



dte -i( u +ir,)t 



2ir ui 



17] 



(21) 



where the limit rj — > has to be taken after the 
integral has been evaluated. 

Thus, the regularization factor ir\ appears ex- 
plicitly in the theory. A possible way to avoid it 
is to use an imaginary-time formulation, which is 
built on a finite contour cut off by a finite temper- 
ature and therefore does not need such a regular- 
ization factor. It is thus one of our goals to clarify 
under what conditions a regularization is not nec- 
essary and justify the use of an imaginary-time 
theory. 



C. Conventional analytic continuation 

In this subsection, we discuss conventional ar- 
guments of the analytic continuation of a real- 
time theory to an imaginary-time theory. We fur- 
thermore illustrate why such deformation of time- 
contour fails for a steady-state nonequilibrium, 
closely following the argument by Doyon and An- 
drei^. 
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In equilibrium, the thermal average of an ob- 
servable A is given as 

with the time-evolution operator S^i,^) = 
e _itjj( tl -t 2 ) wit]l tlie full jjamiltonian # and the 



non-interacting density matrix /5q 



We 



consider that the limit T — > oo exists as discussed 
in the previous section. In the interaction picture 
with Vi(t) = e ltH °Ve~ %tH ° , the above relation can 
be rewritten as 

{A)= ^ Tr Sl (0,-T)p oSl (-T,0) > (23) 



with 



Si(h,t 2 ) = Texp 



dsVi{s) 



(24) 



with the time-ordering operator T defined as the 
time moving in the direction t% — > ti- Using the 
relation, 



—icH\ 



»5/(i) + c,a + c)e lcflo , (25) 



Si(b, a) = e 
one can write 

ft(0, —T)po = poSi(-ip, -ip - T), 



(26) 



in the similar manner as Ref. [2 
written as 



Then (A) is 
TrpoSi(-ip, ~iP - T)S/(-T, 0)4 



^ t" 1 ^ Trp 5/(-i/3,-i/3-T)5/(-T,0) 



lim 



(S I (-ip,-ip-T)S I (-T,0)A) 
(Si(-i0,-ip-T)Si(-T,O)) o 



If we can insert the factor Si(—i(3 — T,—T) 
[denoted as dashed line in Fig. QJa)] between 
Si(—ip, —i(3 — T) and Si{—T, 0), one can close the 
time-contour and analytically continue to the con- 
tour along the imaginary-time (0, —iP) [Fig.QJb)]. 

Using the Wick's theorem and the linked- 
cluster theorem, the terms contributing to (A) are 
of the type 

0, connected j 

(28) 

where the time s = and the interaction times 
{si, • • • , s n } are all interconnected by Wick's con- 
tractions. When the interaction V and the ob- 
servable A are operators local to the QD, one 
can use the relation Eq. ([T^| . We consider a case 



that one of s k in (V/(si) • • • Vj(s n )A)o iCoa belongs 
in the interval [— T, —i/3 — T]. In its connected 
Wick's contractions the operators in A may be 
eventually linked to s& via a forward sequence 
{s' = 0, • ■ • , s' p _ 1 , s' p = Sk} and a backward se- 
quence {s'q = Sk, - " j s 'q-i: s q = 0}- F° r the for- 
ward sequence {s' = 0, • • • , sl_i} with the times 
on the real-axis, we can use Eq. (JTF 



-rmaxdsil,— ,|Sp_ 1 |} < ^9) 



Similar expression holds for the backward se- 
quence. For the last term involving Sk € 
[— T, — i/3 — T], we have a contraction of 
(d(s' 1 ')dt( Sfc ))(d( Sfe )dt( S ;_ 1 )). For -p < 
Im(sfc) < 0, the two factors remain fi- 
nite and give a contribution proportional to 
e-iXlT-M^d+IT-Kl). Therefore, when one 

of the interaction events occurs on the contour 
[— T, — ip — T], the corresponding term becomes 
exponentially small. When traced with local op- 
erator A, the factorization property^ holds 

Sji-ip, -ip - T)Sr(-T, 0) ->■ S T (-ip, 0). (30) 

This shows that the Wick rotation between real- 
time and imaginary-time theory is valid in equi- 
librium and 



(A) 



(S 1 (-tP,0)A)o 
(Si(-iP,0)) 



(31) 



Next we ask whether the same argument can be 
extended to the steady-state noncquilibrium with 
the initial density matrix at time t = —T given by 
po = e -0(#o-*io). i n or der to move p in Eq. (|2"2"1) 
to the leftmost position in the trace, we write H = 
H^ + V* with^ = H -$Y and V"* = V+<f>Y . 
Defining V®(t) = e ltH ° V^e~ ltHo , we can utilize 
(27) the same argument as before to write 



- = (Sn-iP-iP-T)Sf(-T,0)A) o 
[ 1 (Sf(-ip,-ip-T)Sf(-T,0)) • 

(32) 

However, unlike in equilibrium, we cannot use 
Eq. (|16l) for a contraction containing V®(s) since 
V® = V + &Yo contains spatially extended opera- 
tors c ak Ca'k'cr' with contributions well away from 



° Ve~ 



sH, 



° + 



the QD. Furthermore, Vj (s) = 
$y"o with a constant of motion Yq with respect 
to Hq, and Vf'(s) would never lead to an ex- 
ponential decay for the interactions occurring on 
the dashed contour in FIG. QJa). This shows 
that a straightforward analytic continuation of the 
noncquilibrium Kcldysh contour to an imaginary- 
time one is not possible. 
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D. Matsubara voltage 

Recently, one of the authors and Heary22 pro- 
posed that, by introducing a Matsubara term to 
the source-drain voltage, one can extend the equi- 
librium formalism such that the perturbation ex- 
pansion of the imaginary-time Green function can 
be mapped to the Keldysh real-time theory. The 
unperturbed Hamiltonian is written as 



K (itp m ) = H + (iip m - $)Y , 



(33) 



with the Matsubara voltage ip m = Airm / j3 with 
integer m. We take the many-body interaction V 
as perturbation. 

The non-interacting Hamiltonian appears in the 
perturbative expansion in two ways: first in the 
thermal factors e~ l3K °, and second in the time- 
evolution e~ rK ° for the imaginary-time variable 
t e [0, P). The main trick of this formalism is 
that in the thermal factor iy> TO -dependence drops 
out as follows. Since [H ,Y ] = 0, e~ 0it ° = 
e -f3(H -$Y ) e -i Vm @Y _ g ince; with reS p ect t0 t h e 

non-interacting scattering state basis, Y is diago- 
nal and has (half)-integer eigenvalues, e" 1 ' 3 '"' 35 ' = 
1, and we have the important identity 



Po 



(34) 



Therefore, the equivalence of the imaginary-time 
and real-time formalism crucially rests on how the 
double analytic continuation iip m — $ — > and 
t — } it is performed. Since the i<y9 m -dcpcndcncc 
in the thermal factor completely drops out, the 
analytic continuation only concerns the time- 
evolution. For r G [0,/3), e- iVmTf ° ^ 1 and in- 
dependence does not drop out. Thus, one could 
argue that as iip m — $ — > and r it arc taken 
in that order, 



-T[H +(i<p m -*)Yo] 



-tH 



-itH 



(35) 



However, as we will point out in detail later, in- 
tegrals over interaction times may create energy 
denominators of the type (K n — K ln )~ x in the 
perturbation expansions, with K n being the n-th 
eigenvalue of Kq. In such cases, the details of 
the path in the complex plane, along which the 



analytic continuation e v = icp m — $ — > ±i0 + is 
taken, become relevant. On the other hand, in 
the real-time theory, the convergence factor irj in 
the energy denominators determines what poles 
should be chosen. 

III. PERTURBATION EXPANSION 

A. Real-time expansion 

In this section, we investigate under what con- 
ditions the role of the regularization factor 77 of the 
time-independent real-time theory becomes unim- 
portant. We assume that a perturbation expan- 
sion of Eq. ([123")) exists. To illustrate the math- 
ematical structure we choose the fifth-order con- 
tribution (as shown in FIG. [2]) and introduce a 
spectral representation with respect to the non- 
interacting scattering state basis. For the partic- 
ular time-ordering considered in FIG. [2](a) , the 
expression reads 



S a = H) 5 Tr / ds 3 / als 2 

Jo Jo JO 

V I (s 3 )V I (s 2 )V I {s 1 )A 

dt 2 / alhV^t^Vi^po 



dsi 



(36) 



Here we use the notation for intermediate times 
such that U are for the forward contour (— 00 — > 
0, upper time contour) and Sj for the backward 
(lower) contour. We redefine the time as t[ = t\, 
t' 2 = t% — 1\, t\ = tj — etc., and the upper part 
of the Keldysh contour becomes 



dU 



t 

,.0 



dt'o 



—00 




dt'n 



— OO 





dtiVKOVKii+t'a) (37) 
d£ x e iH °^Ve- 



For a spectral representation with respect to en- 
ergy eigenstates, we introduce the convergence 
factor e n ^ tl+t ^ for the reasons discussed in sec- 
tion IA. Then with respect to the non-interacting 
Fock basis |n) and \p), 



J 



r f0 

(-i) a <p| / dt 2 dhVjit^V^n) = 



Vpqyqn 



(E n -E p + ir))(E n -E q + iri) 



(38) 
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A) t 2 U 

O * * * V : O^* X-V, , * ' * O 

X <-* *— * 4 X < * A ~ X *-> X 

53 («0 ^ 31 (b) (c) 

FIG. 2: (a) Keldysh contour in forward direction. Crosses mark interaction points V and the dot an observable 
A. (b) Reversed series of scattering points, (c) Backward Keldysh contour with scattering events equivalent to 

(a) if A is written in terms of QD operators. 



One can do the same for the lower part of the Keldysh contour, 
H)>| / ds 3 / ds 2 / ds 1 V I (s 3 )V I (s 2 )V I (s 1 )\l) 



k (E n -E m - irj)(E n -E k - ir))(E n - E t - irj) ' [ ' 
Therefore the above expression S a can be written as 

g VnmVmkVkl ^ VpqVqw /^qn 

°~ m (E n -E m -iri)(E n -E k -ir))(E n -Et-iri) lp (E n - E p + irf)(E„ - E q + i-q) 1 ' 

Note that all energy denominators consist of one energy anchored at \n) where po acts at t = — oo and 
the other energy of intermediate states \m,k,l,p,q). For the forward contour, the state \n) contributes 
the energy E n + irj in the energy denominator, and E n — ir\ for the backward contour. 

We now consider a counter-time-ordcring as depicted in FIG. [2jb) where the number of scattering 
events on the lower and upper branches are swapped. After an explicit calculation by applying the 
same rules as before, one gets 

rt V n qVqp VlkVkmVmn ( Al) 

6 ~ 4r (E n - E q - ir)){E n - E p - i v ) pl (E n - E x + ir,)(E n -E k + i V ){E n -E m + ir,) Pn ' ( ' 



Starting with the state |n), the numerator 
V nq V qp A p iVikVkmV mn p n in Eq. (gl]) represents 
the reversed process of p n V nm V m kVkiAi p V pq V qn in 
Eq. (gOl). The factor p n V nm V mk V k iA lp V pq V qn is 
understood as the amplitude of the following pro- 
cess 

S a : \n)^\q)^\p)A 

\l) ^ \k) ^ \m) A |n). (42) 

The many-body interaction can be written in 
terms of four scattering state operators as V = 
X) v i2Mi } \' l i } 2^ 1 ^^- With the on-site Coulomb in- 
teraction, 

V = U titzhUg^glgi^l^^l^Ai, (43) 

{a,k} 

where the shorthand notations U = t aj /\/£l, gi = 
gd(ki) and i\)\ a = ij} a . k . a have been used. Note 



that any creation of a particle ipj is associated 
with the factor tig* , and the annihilation iftj with 
tjgj. For the observable A we consider a one-body 
operator A = ^ ai2'0iV'2 for simplicity. The oper- 
ator V creates up to two particle-hole pairs of type 
V>, and for a non-zero matrix element (n|V|m), 
\n) and \m) differ only by up to one particle-hole 
pair per spin channel. Thus, in the above pro- 
cess Eq. (|4"2"j) . which starts and ends with |n), the 
product of creation operators </>„ fccr must match 
the that of annihilation operators ip a ka- There- 
fore, the matrix element for the process Eq. (|4"2"j) 
must be of the form 

S a : \tigi\ 2 \t 2 g2\ 2 ■ ■ -tigidijtjgj. (44) 

Similarly, the process for Sb-term 

S b : \n) — » | to) — > \k) — > 

10 A \p) A \q) A \n) (45) 
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must contain the same set of {ip^,ip} with the 
same states, only in the reversed order. The ma- 
trix element for the process then becomes 

Sb ■ I 2 1*2.92 1 2 ■ ■ •'.«.';,'',•''.'/,'• (46) 

If the operator A satisfies the following property 

9d(ki)aij[g d (kj)]* = g d (k ] )a jl [g d (k l )]* , (47) 

the matrix elements for counter-contours (a) and 
(b) match, i.e. 

VnmVmkVklAipVpqVgn = V nq V qp A p iVlkVkmVmn- 

(48) 

With this condition, S a (rj) = Sb{— rj), and S a + Sb, 
inside the expression for (A) , is independent of the 
sign of rj and has a well-defined limit of ?y — > ±0. 
The above argument can be repeated for any order 
of the perturbation expansion, i.e. the use of a 
spectral representation is permitted and the result 
independent of the convergence factor rj provided 
that the contour has itself as the counter-contour, 

S a (v) = Sa(-v)- 

Which of the physically interesting operators do 
satisfy the above condition Eq. (|4"T)) respectively 
(|48l) ? It is easy to see that it is true for any oper- 
ator A which is a simple function of n da — d\d a . 
A general two-body operator 

A = y] 01234^1^2^3^4 
1234 

also falls into this class if it satisfies 



gd(h)gd(kj)aij 

III!) [gd(k n )gd(k m )}* 
9d{ k n )gd{ k m ) &nmij [9d(ki)g d (kj)]* ■ 



(49) 



Unfortunately, the current operator Eq. (fT4"f does 
not satisfy the condition Eq. (|4"T]) . and a direct 
analytic continuation is not available, as we will 
discuss shortly. Therefore, we have to resort to 
the Mcir-Wingreen formula^ which relates the 
current to the spectral function. 

We have so far ignored coinciding energy de- 
nominators in the perturbation expansion lead- 
ing to overlapping (5-functions. For the sake of 
simplicity we consider a second-order contribu- 
tion from Eq. (|23p . By expanding it into different 
time-orderings, we obtain 

dh [ dt 2 p Vi(t 2 )Vi(h)A 

t-T 

dh / dt 2 Vi(ti)p Vi(t 2 )A 



(50) 



T JO 

1-0 pt! 

' dh / dt 2 Vi(h)Vi{t 2 )pQA. 

T JT 



We now introduce the convergence factor e vt and 
take T — > oo to obtain the expression 



E 

nml 



Pn 



(E n - E m + ir])(E n - Et + ir]) 

Pin 

(E m - E n + ir))(E m -Ei- irj) 
Pi 



{Ei- E n -n 1 ){E l - E m -i V ) 



VnmVmlAl n . 



which needs precaution when the two energies 
in the denominators become equal, because the 
contribution will be a product of two <5-functions 
with the same argument. One must be careful 
when one performs the limit T — > oo. To see 
this let us go back to the timc-dependent descrip- 
tion. By keeping T finite, contributions of the 
form S(E n — E m ) 2 will actually amount to terms 
proportional to T 2 from the integrals. Combining 
all three integrals we obtain the coefficient to the 
T 2 -term (i.e. <5 2 -term) proportional to 



(pn - 2Pm + Pl)V nm V m .iAi n X 



nml 



6(E n - Em)6{E m - El) 



(51) 



In equilibrium p n = p m = pi for E n = E m = E t 
and this term vanishes identically. The argument 
can be easily extended to arbitrary orders in the 
perturbation expansion. 

In the case of noncquilibrium the situation is 
more complex. Here we discuss in detail what 
happens to Eq. (fSTj) . Wc consider the case \n) ^ 
| to) 7^ |/), while E n — E m — Ei. Suppressing the 
(5-functions, Eq. ((BTj) has the form 



e -^n (e 0«W _ 2e ^o m + e ^)V nmVm lA ln . 

In the matrix element V nm VmiAin, the transition 
\n) — > \m) — > \l) — > \n) involves a certain series of 
particle-hole excitations. For instance, \n) —>■ \m) 
is given by an exchange of two particlc-holc pairs, 
^a 1 k 1 J- , a2k 2 a^ a3k3a ,^ ai k i a' in V, and similarly 
for |m) — > |/) and \l) — > \n). However, since any 

creation of i^ aka should be matched by ip a k<j only 
up to 6 indices are independent. Given a par- 
ticular set of the 6 indices of wave-vectors and 
spins {fci<7i, k 2 a 2 , ■ ■ ■ ,fcg(76}, different permuta- 
tions of the above 6 pairs of {V'fc o-- > V'fcio-i} m 
VVA determines the matrix element V nm VmiAi n . 
Now, we sum over all possible combinations of 
reservoir indices {«!,•■■ ,ag} (while keeping the 
k- indices unchanged) for the all twelve {ifj,ip} 
operators. The matrix element V nm VmiAi n oc 
ll l =i,6*aj5(efc I )| 2 - Since the product of \g(e kz )\ 2 
are invariant, we collect all possible reservoir 



weights in Y\ i=1 6 t\ e^**^™." 1 .'} and each of the 
three sums in Eq. (|5ip become the same, i.e. the 
whole contribution vanishes. A detailed discus- 
sion of the mathematics can be found in Ap- 
pendix [A] 

In summary, if the observable A satisfies 
Eq. (|47|). the energy integration in the pertur- 
bation expansion can be interpreted as principal- 
valued, similarly to equilibrium. In Appendix [B] 
we provide as an example the fourth-order contri- 
bution to the QD-electrons self-energy and show 
explicitly that the above properties are satisfied. 
Since the structures appearing in higher order are 
of the same type as discussed above, we may ac- 
tually infer that this property holds in any order 
of the perturbation expansion. 



This expression has the same mathematical struc- 
ture as in the real-time theory. Even though 
we considered only one time-ordering in the 
imaginary-time theory, the upper and lower in- 
tegral limits in dr dr' combine to create the 
same permutation of terms as in the real-time the- 
ory^. 

We have seen earlier that, in the real-time the- 
ory, energy denominators can be interpreted as 
principal- valued since all (5-function contributions 
from the energy poles vanish. Therefore, if we 
interpret the energy denominators as principal- 
valued as i(p m —> <f> 



B. Imaginary-time expansion 



K n -K n 



V 



E n — E„ 



(55) 



Unlike the real-time theory, the imaginary-time 
description is formulated on a finite time interval 
of [0, /3), and there is no need for a convergence 
factor e 1 '*. Therefore, the energy integrals appear- 
ing in the equilibrium theory are always principal- 
value integrals, which we confirmed in the previ- 
ous section HTCl 

In nonequilibrium, with the imaginary-time ef- 
fective Hamiltonian K(itp rn ) = Hq + e^Yo + V 
$), the thermal average is defined as 



{A) 



Tve-^A 



Tj- e -0K 

The Boltzmann factor can be expanded as 



-PK 



T T exp 



P 



drVf(r) 



(52) 



(53) 



with Vi(t) = e TKo Ve- TKo V and T T denoting the 
time-ordering operator for re [0 — > f3]. We con- 
sider a second order expansion to understand its 
mathematical structure, 

Tre -Mo f f dT'V I {T)V I {T')A 



dr / dr' 
Jo 



V« e r(K n -K m )y T'(K m -K x ) v 



E 



(K n -K m ){K n -Ki) 

Prn 

'(K m -K l )(K m -K n ) 
Pi 



the terms in the imaginary-time theory indeed 
match those of the real-time approach. 



In section IIV A 11 we calculate the double oc- 
cupancy from continuous-time quantum Monte 
Carlo method, and numerically verify that the 
analytic continuation procedure outlined so far 
works accurately in all orders of perturbation 
theory as well as for the resummcd perturbation 



C. Single-particle self-energy 



The analytic properties discussed so far can be 
used to examine the single-particle self-energy for 
the Anderson impurity model. The imaginary- 
time second-order self-energy in the Coulomb in- 
teraction U can be written as22 



Z (2 \iLj n ,e v ) = V / de- 
^ J i 



cr 7 (e) 



(56) 



(Ki - K n )(Ki - K„ 



VnmVmiAin- (54) with the spectral function 



10 



CT 7 (CJ) = U 2 



[ / deiA (ei) 



- f*)h + (1 - /i)/a(l - /a)] S(u - a + e 2 - ea) (57) 



a 1 +a 2 +"3=7 



for the 7-branch cut (7 = ±1,±3), where 

. , . r/7r 

(e - e a ) z + T- 

dcnotes the non-interacting spectral function of 
the QD level and f a = [1 + e -^(«-«*/2)]-i the 
Fermi-Dirac factor for the a-th reservoir. 

Recently, it has been proposed^ that an in- 
clusion of higher-order contributions will mainly 
modify the spectral function cr 7 (e), leading to a 
e v dependence like 

£(io,„, e v ) = J2 [ de . ■ ( 58 ) 

Based on this expression, one can try to fit 
<r 7 (e, €(p) to the numerical single-particle self- 
energy generated from quantum Monte Carlo cal- 
culations. However, in order to establish the 
existence of an analytic continuation limit of the 
imaginary-time self-energy, one should first show 
that the real-time self-energy possesses the ana- 
lytic property discussed in the previous section, 
namely that the energy poles are principal- valued. 
The rather lengthy and technical argument is pro- 
vided in Appendix [B] for the fourth-order self- 
energy diagrams. It can be shown explicitly that 
contributions involving products of (5-functions 
with identical argument vanish identically, result- 
ing in the necessary analytic properties discussed 
in the previous section. 

Again, investigating the general structures ap- 
pearing in the perturbation expansion of the self- 
energy, we are confident that this property indeed 
holds in any order and also survives the resum- 
mation of the series. The latter aspect, however, 
cannot be proven rigorously, but is strongly sup- 
ported by the numerical evidence from our Monte- 
Carlo simulations. 

In a recent work by Dirks et al^ and a ac- 
companying paper to this work, a general analytic 
continuation approach based on the multi- variable 
complex function theory and its double analytic 
continuation of (iu> ni iip m ) have been systemati- 
cally studied. 



D. Forward and backward steady-state 

We have seen in Section IIII Al that we need 
Eq. (|48p for any sequence of matrix elements in 
order to establish the equivalence of the real and 
imaginary-time theory. In order to close the for- 
mal discussions, let us re-examine the complex 
conjugate of the matrix elements in relation to 
the forward- and backward-in-time propagation of 
scattering state density matrix. 

Assume that we propagate a non-interacting 
density matrix po = cxp[— /3(Hq — &Yq)] from the 
initial time t = —T to the present in the for- 
ward direction. Then, according to Gell-Mann 
and Goldberger— , we obtain 

Pout = V r 'e' lCT (e^Tp o)e -nT dT 

= 77 / e^ CT p e^ T dT 
Jo 

V - 

= — r~rP° 

77 + iL 

= Po H -^—C v po , (59) 

-C + ir\ 

with Cv the Liouvillian representing the interac- 
tion parts not contained in Cq. p ou t is the fully 
interacting density matrix at t = and po non- 
interacting density matrix at t = 0. The mean- 
ing of the above equation is that we unwind a 
non-interacting density matrix to a remote time 
t = —T and re-evolve it with full interaction to 
the present time. By taking the average over the 
remote time T, we filter out transient oscillations. 

Alternatively, we can also consider a backward 
propagation of density matrix evolving from the 
remote future by writing 

Pin = V / e* CT (e-^> T p )e^ T dT 
Jo 

= Po H -p. — —Cvpo- (60) 

— L — irj 

If we initially choose po as the density matrix 
of a quantum dot system of disconnected dot and 
reservoirs, £y = £ t + Cjj includes both the hop- 
ping to the leads and the Coulomb interaction on 
the dot. We first construct the scattering states 
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with respect to the hopping, and then with re- 
spect to the Coulomb interaction. After the first 
step, the scattering states become^! 

i>lk„,out = C aka + jfi9d(k)4 +- (61) 

^U<» = 4, + ^(^4 + ..-,(62) 

and we can construct respective scattering-state 
density matrices pot,out and pot, in with Cy — Cjj- 
The coefficients appearing in front of the dot op- 



erators d\,d a etc. for the out and in-scattering 
states are the complex conjugate of each other. 
Therefore, the matrix elements of the interaction 
V = Un dt n di , written in terms of Vw,{ottt,m}" 
basis, are complex conjugate to each other, i.e. 
V nm = Vim (with the tilde denoting the in- 
scattering basis). 

We can now repeat the arguments from Sec- 
tion IIII Al for the backward propagation of the 
density matrix as shown in FIG. [U[c) and find 



c £—< (PL - R_ 4- - R_ 4- i<n\ & 



(E n -E q + irf){E n -E p + i V ) p 1 {E n - E t - ir))(E n - E k - ir))(E n - E m - in) rn ' 

For observables satisfying A nm = A*^, this expression becomes identical to S a in Eq. (|4"0j) . The 
same argument holds in any order of the perturbation expansion, and we have Tr Ap out = Tr Ap~i n and 
(A) = \ {{A) out + (A) i n ). Therefore, from Eqs. ([55]). 0, we have 

(A) = (i) + (ii (-^ + Cypo) = (i)o + (AV ( JL) Cypo) , (63) 



i.e., the conditions for replacing the energy de- 
nominators by their principal- values, as discussed 
in section llll A[ correspond to a measurement pro- 
tocol where the observable A has the same ex- 
pectation values with respect to the forward- and 
backward-propagating density matrices. 



IV. STATIC EXPECTATION VALUES 
A. Theoretical background 

We have shown that steady-state expectation 
values of certain local observables A can be ob- 
tained from analytical continuation of expecta- 
tion values calculated within the imaginary time 
Matsubara-voltage formalism. As long as we 
know the analytic structure of these objects, this 
can be done easily However, for a model with true 
two-particle interactions, one eventually has to re- 
sort to numerical evaluations, and an analytical 
continuation in general requires a more involved 
computational technique. We therefore want to 
provide in the following a representation which 
allows the use of standard tools from equilibrium 
many-body theory. 

A numerical method gives (A)(iip rn ) and let 
(A)(z lp ) be its analytic continuation. We may 



write formally 

(A){Z V ) = (i)const + Xa{z v ) (64) 

where the part Xa{z) is holomorphic in the upper 
and lower half plane, with singularities only on 
the real axis. If one can furthermore show that 
the zxa(z) is non-singular in the limit z v — > oo, 
one can finally infer that a spectral representation 
with respect to the jump function on the real axis 
exists and hence 

(A)(i<p m ) = (i)const+ / T. gAi l\ d<p 

J [l(f m - <P) - Lf 

(65) 

Note that the latter property is not necessarily 
guaranteed and has to be proven individually for 
each observable. 

Once the validity of the representation (|65[) is 
established, one only needs to obtain the "spec- 
tral function" qa (y). One evident method to 
calculate the Matsubara voltage data (A)(itp m ) 
for the observable A with respect to the effective 
system with non-hermitian Hamiltonian at Mat- 
subara voltage i(f m is via a QMC simulation.— 
For such data with statistical noise, one then 
typically employs a maximum-entropy approach 
(MaxEnt)4£ The implementation of a MaxEnt 
estimator for the physical expectation value is 
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rather straightforward. The values for differ- 
ent itp m are truly statistically independent, and 
only the variance and correlation between imagi- 
nary and real parts of a single iip m value play a 
role. However, one still needs accurate and un- 
biased measurements of imaginary-voltage data 
over a large range of ip m ^ This latter require- 
ment makes the use of a continuous-time quan- 
tum Monte-Carlo (CT-QMC) algorithm manda- 
tory. In particular, the necessary estimation of 
the constant offset (A) const in Eq. (|65|) is possible 
only with CT-QMC, because at present no direct 
measurement algorithm for this quantity is avail- 
able and one must determine it from the tail of 
(A)(iip m ) by fitting it to 



{A){lip m ) {A} const + 



CA 



{ifm)'' 



+ 



(66) 

In practice, a weighted least-square fit yields reli- 
able values and error bars for (^4) const- Via Gaus- 
sian error propagation it is then possible to incor- 
porate the uncertainty of (A) const into the covari- 
ance matrix of the quantity (A)(iip, n )~ (A) const £2- 
In general, the spectral function qa{^p) needs 
not to be positive semidefinite, or show any sym- 
metry relations with respect to ip. Since on the 
other hand the MaxEnt method is only applica- 
ble for the inference of positive definite functions, 
a shift function g s hi{t( l p) of the spectral function 
Qa(}P) has to be introduced, which makes the to- 
be-inferred q'a(<p) = qa^) — Qshnt(p) positive. We 
also employ a symmetry condition 



Pshift^) = £>shift (-(f), 



(67) 



because this choice is robust with respect to the 
physical result 



Qa(<p) 



a = ±l 
(A) const 



(68) 



dip ■ 



In the following we want to prove that the dou- 
ble occupancy or magnetization obey this con- 
straint, i.e. have a representation, where (A) const 
is a real number, and qa{ip) € R is a real- valued 
spectral function. 



1. Double Occupancy 

The double occupancy in Matsubara-voltage 
representation is defined as 



D(iip m ) := (n d ^n dti ) 



K(iip„ 



(69) 



where the expectation value is taken with respect 
to the m-th effective equilibrium system. 

We will first show that the representation (|6"5j) 
holds for the double occupancy, i.e. that we have 
indeed 



D(itp m ) = D + / dip - 



tip 7] 



$ - tp 



(70) 



We restrict the discussion to the case of particle- 
hole symmetry and symmetric coupling to the 
leads, Tl = Tr. Within the Matsubara-voltage 
approach, one can - for fixed iip m - employ the 
standard techniques of equilibrium many-body 
theory and obtains the standard result^ 

D(iip m ) =(n f )(ni) 

+ «77 X! iu n )G(i(p m ; iw n )e lu)nV . 

t^n 

(71) 

Due to particle-hole symmetry, we have 
(n^)(ni) = 1/4. Furthermore, from the dis- 
cussion in section IIII CI we can infer that at 
least the Green's function decays like \/iip m 
and hence allows for the existence of a spectral 
representation (|70|) . as long as there is only a 



single branch cut at Im z v = 0. 

The real-valuedness of spectral function and 
constant offset remain to be shown. The general 
relation G(~itp m , —iuj n )* = G(itp m , iu> n ) holds for 
Green's function and self-energy Inserting this 
into Eq. (JTTJ) , we find 



D(-icp m )* = D(iip m ). 



(72) 



Consequently, the real part of D(iip m ) — D{— iip m ) 
vanishes. Using the symmetric coupling to the 
leads, we have an invariance of the Green's func- 
tion and self-energy under {iip m — $) <-> — {itp m — 
<f>). As a result, Do is an actual constant which 
is obtained for both, upper and lower half plane. 
Due to the symmetry of Im D(itp m ), Dq is real. 
By inserting the representation ([TO"]) into Eq. (jT2"j) 
we also see that £>d(w) is real- valued. 

For example, let us consider the equilibrium 
setup, i.e. $ = 0. At half filling and symmetric 
coupling to the leads, the function 

Re£>$ =0 (i</'m) = Re-D^oMVm), (73) 
Im£»* =0 (^ m ) ee 0. (74) 

This is compatible with a conventional bosonic 
spectral representation 



D<f> =0 (itp m ) 



dtp-?^ + Do, (75) 

iip m - tp 
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with an antisymmetric spectral function 

Qd(<P) = -Qd{-<p)\ to(^>0)<0 (76) 

and the offset Do > 0. Eq. (|74|) is not evident for 
asymmetric couplings or off particle-hole symme- 
try, because here Go(itp m ,iT) is not real. 



2. Magnetic Susceptibility 

An observable which is much more sensitive 
to the Kondo effect is the magnetization M := 
((n^) — (nj.)) in the presence of a magnetic field B 
in z-dircction respectively the magnetic suscep- 
tibility x = M/B of the quantum dot, because 
it directly probes the spin degree of freedom of 
the dot electrons. In equilibrium, a strong depen- 
dence on the temperature is observed, on the scale 
of the Kondo temperature.— 

As for the double occupancy, the validity of a 
spectral representation 

M{i Vm ) = M + [ d^ . gM ^ } (77) 

J vp m - <P - ip 

can readily be confirmed. Starting from the sym- 
metry G(—iip m ,—iu> n )* = G(itp m ,iu) n ), one can 
again show that M(—i<p m )* — M(iip m ), and the 
same arguments apply concerning the interchange 
{iip m - $) O -(i(p m - $)• 



B. Numerical effective-equilibrium data 

Let us now turn to the discussion of actual 
numerical data for magnetization and double oc- 
cupancy from the quantum Monte-Carlo simula- 
tions. As the first step, we analyze these data 
with respect to the auxiliary variable tp m , and 
want to argue that they have a physical interpre- 
tation with respect to the actual voltage <£>. In 
particular, the convergence of the numerical pro- 
cedures described below implies full consistency 
of the Matsubara-voltage formalism with regard 
to the numerical data. 

We find that effective-equilibrium data come 
along with characteristic energy scales which - af- 
ter analytic continuation - may translate almost 
directly into energy scales with respect to the ac- 
tual source-drain voltage $. It is therefore worth- 
while to discuss the dependence of the effective- 
equilibrium expectation values as a function of ip m 
for given physical parameters /?, U, and $. 

a. Dependence on <f>. The first thing to no- 
tice is that the dependence of the shape of the 



curves M(iip m ) and D(i<p m ) on $ is rather mod- 
erate: for the examples considered, we do not ob- 
serve any new characteristic energy scales with 
respect to the Matsubara voltage <p m emerging or 
disappearing as a function of the physical voltage 
$. The most striking influence of $ is a change 
of the offset of the curves Do and Mo- The offset 
is changed monotonically function of $ and 
cannot explain features such as dips and peaks 
which are found in the analytically continued data 
(cf. next section). This is the very reason of our 
claim that low- to intermediate-energy scales with 
respect to ip m rather directly translate into low- to 
intermediate-energy scales with respect to 3>, al- 
though tp m has no direct physical meaning itself. 

Let us substantiate the above statement by the 
data plotted in Figs. [3l and l4al In Fig.[3j effective- 
equilibrium double occupancy curves are shown 
over a wide range of values of the physical voltage 
and Coulomb interaction. Each curve exhibits a 
dip at ip rn = 0. As already pointed out above, 
the dependence on $ is rather mild, except for 
the offset. The same behavior is observed for 
the magnetization in Fig. 0aj i.e. the voltage <f> 
merely introduces an overall shift and a moderate 
smoothening of the structures. 

b. Limiting behaviour tp m — > ±oo. For each 
U and $ a different limit Do is obtained as (p m —> 
oo. If the values /3, U, and in particular ip m 
arc large, the effective-equilibrium QMC simula- 
tions start to suffer from a significant sign prob- 
lem. This may result in particularly noisy tails 
such as the ones for the data with largest $ in 
figure I3dl In these cases, the estimate of Do is 
subject to much uncertainty and limits the statis- 
tical accuracy of physical expectation values. 

c. Dependence on U . As U is increased, the 
depth of the dips in the double occupancy curves 
also increases. On the other hand, neither the 
width nor the shape change significantly. In par- 
ticular, the emergence of a Kondo scale Tk can- 
not be inferred from these data. Interestingly, for 
small U, the relative contribution of the constant 
term Do is large compared to the height of the 
peak which emerges around tp m ~ 0. As the inter- 
action increases, the central peak becomes more 
pronounced, and the physical expectation value 
increasingly depends on the structure of the peak. 

For the magnetization in Fig. I4bl a similar 
picture seems to emerge at first glance, namely 
a strong increase of the offset Mo with U to- 
gether with a more pronounced peak structure at 
ip m = 0. The strong increase of both is readily 
understood as with increasing U the system forms 
a local moment which is aligns with the external 
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FIG. 3: (color online) Real part of the effective-equilibrium double occupancy as a function of the 
Matsubara voltage ip m at several values of interaction strength U and bias voltage <£>. 



field. 

d. Kondo effect. Up to now there seems to 
be no evidence whatsoever for the presence of the 
Kondo scale 7k in the data presented so far. On 
the other hand, the generation of this many-body 
scale is usually considered as crucial test for any 
method proposed for studying the Anderson im- 
purity model. As already pointed out, it is quite 
apparent from the data in Figs. [3] that 7k ob- 
viously does not appear to be relevant for this 
quantity; a fact that is already well known in equi- 
librium. There the scale 7k shows up only in a 
very indirect way as renormalization of the zero 
temperature value respectively the scale regulat- 
ing the approach to \t£^ 

The situation is different for the magnetization. 
Here, the Kondo scale plays a crucial role^ as it 
determines the field-strength necessary to break 



up the Kondo singlet. Hence it must show up in 
the magnetization; in particular, one must actu- 
ally expect a scaling behavior with Tk for small 
enough fields. Let us therefore plot the mag- 
netization as function of Matsubara voltage in 
the form M(<^ m /7k) for values of U beyond the 
weak-coupling regime for fields and voltages much 
smaller that the corresponding equilibrium Kondo 
scales. The result is shown in Fig. [5j Evidently, 
the width of the peak in the effective-equilibrium 
magnetization data is nicely scaling with the equi- 
librium Kondo temperature, i.e. for different val- 
ues of U the peak structure is essentially left in- 
variant at fixed values of B, <f> and T. 
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C. Results for real voltages 

In this section we will introduce the MaxEnt 
procedure used to infer the spectral functions 
qd{<p) an d Pm(v 3 ) from the effective-equilibrium 
QMC data. Based on this analytical continuation, 
we then will discuss the physical results obtained 
from the auxiliary Matsubara voltage data. 



1. MaxEnt procedure 

Based on the effective-equilibrium data and the 
exact relation (|65|) . it is in principle possible to 
uniquely reconstruct the spectral function qa{<p) 
and the offset (A) const . This is almost completely 
analogous to the conventional Wick rotation. 



However, because in practice a finite set of data 
is considered, the inversion of equation (|65[) is 
no longer unique. On top of this, the quantum 
Monte-Carlo data are not exact but merely Gaus- 
sian random variables. One may easily verify that 
the noise associated to the variables is amplified 
by the inversion of equation (p5|) . As a conse- 
quence, it will always be possible to find quali- 
tatively very different functions qa^) which are 
in agreement with the QMC data. In particu- 
lar, these functions will yield physically different 
predictions via equation (|68[) . The problem to ob- 
tain physical results from the effective-equilibrium 
data is thus ill-posed. 

Since essentially the same integral equation (|65[) 
also relates imaginary-time and real-time proper- 
ties of conventional Green's functions, this issue is 
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well-known to the community4£ Although no so- 
lution to the problem can be provided, Bayesian 
inference provides a framework to systematically 
incorporate a-priori information about a quantity 
into an estimate. The estimate is most likely 
with regard to the prior information at hand. 
The resulting method is called Maximum Entropy 
(MaxEnt) 4£ 

Let us consider the situation in which the offset 
(^4)const has already been determined via a least- 
square fit. Via error propagation it has been pos- 
sible to determine the covariance matrix of the 
quantity (A) — (A) const , i.e. the imaginary- voltage 
values of the quantity xa(z v ) in equation (|G4l) . 
The remaining task of the MaxEnt is to infer the 
spectral function g^ip). Let us furthermore as- 
sume that the data have been sufficiently trans- 
formed with a shift function, such that the func- 
tion 



(78) 



is positive (see section HV A[) . 

The default model for g' A (<p) is then a posi- 
tive definite function which in principle should 
contain features which determine in particular 
the high-energy behaviour, if known4£ In the 
case of Green's functions, perturbation theory or 
higher-temperature solutions often give good de- 
fault models4£ In our case, apart from that we 
used a shift function to construct the positive 
spectrum, nothing is known about the function, 
so a flat default model is preferable. As conse- 
quence, we use the shift function itself as the de- 
fault model in the actual computation. For sim- 
plicity, let us call the to-be-inferred spectrum g(y>) 
and the default model Pdof(y)- 

On the one hand, the default model gives rise 
to a relative entropy^ 



S 



dip 



q(<p) - Qdef(f) " Q{f) log 



gOg) 

£>dcf(v?) 



of the spectral function. On the other hand, 
the (transformed) effective-equilibrium simulation 
data with mean values 5j and covariance CV,- yield 
the measure 



^ Nqmc 



id 



{a, 



Vi)- (79) 



for the quality of the fit. Here yi are the fit val- 
ues which result from transforming the considered 
g{<p) to the data space, and Nqmc is the number 
of QMC data points a,. Within the MaxEnt it fol- 
lows that a functional Q — \ 2 — ceS must be min- 
imized, where a > is some hyperparameter4£ 



In order to determine a, there are several meth- 
ods, for example the "historic" and the "classic" 
MaxEnt M> The former extracts information from 
the Monte-Carlo data up to the point at which the 
X 2 = Nqmc, Le. the MaxEnt rcgularization pa- 
rameter is fixed to the value at which x 2 = Nqmc- 
The latter ("classic" MaxEnt) extracts informa- 
tion from QMC data to a larger extent. Based 
on the probability distribution implied by the de- 
fault model and maximum-likelihood functionals, 
a posterior probability of the MaxEnt rcgulariza- 
tion parameter a is maximized. Because informa- 
tion from the default model is again incorporated 
rather explicitly, this strategy is particularly good 
for default models which are close to the actual so- 
lution. A rather general feature of "classic" Max- 
Ent appears to be that the x 2 value of the inferred 
estimate is generally much smaller than the "his- 
toric" value of Nqmc- Our feeling is that this 
aspect makes the "classic" estimate more sensi- 
tive to statistical fluctuations and vulnerable for 
over-fitting, but on the same side, the estimate 
is less biased. A similar increase in fluctuations 
was pointed out in a recent study^ At least if 
Bayesian evidence coming from the data is weak, 
the "historic" MaxEnt, on the other hand is more 
biased towards the default model value, since its 
estimate is more conservative with regard to the 
X 2 ■ In our case, the default-model estimate is 
given by the constant offset Do, because our de- 
fault models are chosen to be even functions with 
respect to ip. 

As shift functions, wide Gaussians with width 
a = T were used, i.e. 



£?shift 



(ip) = X- e~*~l 2a2 



(80) 



The amplitude of the functions was varied in such 
a way that positive functions could be inferred. 
The different values for differently scaled func- 
tions give rise to a certain interval of expectation 
values, which will be plotted result, in the 
following. An example for the set of inferred func- 
tions obtained for a single non-equilibrium system 
is shown in figure HI The left panel shows the ac- 
tually performed MaxEnt for the shifted spectral 
functions, using "historic" MaxEnt. Resulting 
from a flat default model for the function qd(!P), 
the shift function acts as default model here. In 
this case, choosing a parameter A < 0.01 yields ar- 
tifacts in the physical solutions, because the neg- 
ative regions of g(tp) cannot be represented any 
more. The corresponding actual spectral func- 
tions g(<p) , obtained by subtracting the shift func- 
tion ([50)) from the data in the left panel, are shown 
in the right panel of Fig. [6] The flat default 
model represents our lack of prior information 
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FIG. 6: (color online) MaxEnt inference process for the double occupancy. Parameters are U = 5, 
e$ = 0.25r, f3 = 20r~ 1 . Due to lack of prior knowledge, we use a flat default model, i.e. the shift 
function p s hift (<^) , see Eq. (pi)) . Remember that the actual spectral function Qd^) was shifted to a 
positive one, g' D (ip), via equation (|78j) . One finds that the different equivalent ways of imposing a 
flat default model for Qd^) yield practically the same spectral function. Nevertheless, computing the 
physical value (|68|) yields values which are distributed over a certain range. This range is displayed as 

error bars in the results plots Figs. [7] and [5] 



about the solution and the preference of a smooth 
solution in case of uncertainty. In general, the dif- 
ferent realizations of a flat default model with the 
shift functions yields almost but not exactly the 
same spectral functions. In case of limited QMC 
data quality, it is well knownii that the usage of 
a flat default model yields less accurate spectra 
than an appropriately constructed more informa- 
tive default model. For example, in case of con- 
ventional equilibrium spectral functions of Fermi 
or Bose systems, a default model should prefer- 
ably obtain the correct low-order moments, which 
can often be computed exactly. It can thus be ex- 
pected that quantities that arc calculated from 
the spectra inferred using the flat default model 
are biased towards a certain value. Nevertheless, 
an increase in data quality will eventually reduce 
the bias of the estimated quantity. We also expect 
that the precision of our method can be increased 
by the development of default models which con- 
tain additional information like moments. How- 
ever, at present such type of information is not 
yet available. 

In order to obtain a rough estimate on the error 
of a physical estimate, we will plot the intervals 
which are generated by computing the estimates 
for different values of A. Typically, a range from 
A = 0.01 to A = 0.16 is imposed, unless the neg- 
ative regions of g(<p) cannot be represented. For 
the magnetic susceptibility, the same strategy is 



used. 



2. Double occupancy 

We will now discuss the analytically continued 
data of the double occupancy and compare it with 
respect to zero-temperature second-order pertur- 
bation theory^ In figure [7] we show double oc- 
cupancy data for different values of the Coulomb 
interaction computed with the two different Max- 
Ent estimators. 

The complementary behaviour of the two esti- 
mators may be well observed in Fig. [7J In the 
large-bias limit, in which the perturbation theory 
may be expected to be correct, the classic estima- 
tor is closer, and the historic estimate is systemat- 
ically too high. This is in agreement with our ex- 
pectation that the historic estimate will be biased 
from above in case of rather weak Bayesian evi- 
dence from QMC data, because the ill-posed con- 
tinuation problem is particularly severe at high 
energies 4i Apart from some fluctuations in the 
"classic" estimator, the same curves are predicted 
for small voltages. It is important to note that 
error bars in the figures do not denote statistical 
errors (which cannot be estimated) , but the range 
of values which a given set of symmetric default 
models generates. 

As compared to the second-order perturbation 
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function of the bias voltage at different values of 
[/, as compared to second-order perturbation the- 
ory. In addition, the dashed lines show the tem- 
perature dependence of (n^ni) in equlibrium as 
obtained by NRG, assuming e$ = k^T (see text). 



theory, we find that both methods agree per- 
fectly for interaction strength U = 3T. Also both 
methods predict a minimum in the double occu- 
pancy at voltage e<E> sa 21? which slowly shifts to 
larger values of $ and becomes increasingly distin- 
guished as the interaction is increased. There is, 
however, a clear difference concerning the magni- 
tude of this minimum, which appears much more 
pronounced in the QMC data as in the perturba- 
tion theory. Note that this seems to be the case 
for both MaxEnt estimators. At present the ori- 
gin of the deviation is not clear. 

One of the issues related to the $ dependence 
of stationary non-equilibrium quantities is to what 
extent they can be mapped onto an effective equi- 
librium temperature dependence. To have an idea 
whether this mapping works, we included in Fig. [7] 
also the corresponding curves for (n-|-n^)(T) as ob- 
tained from an NRG equilibrium calculation, as- 
suming e$ = k^T. Quite apparently, the val- 
ues at $ — > nicely coincide, which also tells us 
that the Matsubara voltage QMC reproduces the 
proper low bias results even for strong coupling. 
Note that perturbation theory here deviates sys- 
tematically with increasing U. However, the de- 
pendence of (n-^n^)($) cannot be mapped even 
qualitatively onto (n^ni) (T) by a simple ansatz 
<f>=a ■ T with some value a for any of the U values 
considered here. From this observation we would 
thus conclude that such a mapping is - at least for 
the simplest possible quantity - not appropriate. 
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FIG. 8: (color online) Magnetic susceptibility as 
a function of bias voltage in the Kondo regime 
U = &T at f i B B = k B T K /2, T = T K /2. The dot- 
dashed line represents an equilibrium NRG calcu- 
lation for T > Tk/2, rescaled in both magnitude 
and temperature to match the low-bias behavior 
of historic MaxEnt (see inset). The double-dot- 
dash curve finally is a fit of historic MaxEnt to 
some scaling function (see text). 



3. Magnetic Susceptibility 

Similarly, the magnetic susceptibility may be 
computed as a function of the bias voltage by 
analytical continuation of the QMC data. As 
an example, we show the result for U = 8r at 
the temperature T = Tk/2 and magnetic field 
I-LbB = UbTk/2 in Fig. [8] When we compare 
our continuation results at $ — >• to the ex- 
act low-bias limit (i.e. the equilibrium value, dis- 
played as a cross in Fig. the historic MaxEnt is 
again more strongly biased than the classic Max- 
Ent, i.e. the deviation from the equilibrium value 
is stronger. With insufficient QMC information, 
the outcome is more biased towards the flat de- 
fault model and from Eq. (|55|) the integral van- 
ishes in such limit. The constant offset Mq lies 
below the actual physical limit, and therefore, as 
QMC quality improves, our estimate approaches 
the correct limit from below. Again, the classic 
MaxEnt is subject to stronger fluctuations. 

In physical terms, the decay in magnetic sus- 
ceptibility is because of the destruction of the 
Kondo effect due to the decoherence introduced 
by the bias voltage. This is in principle sim- 
ilar to the equilibrium behaviour found as a 
function of temperature^ The scale on which 
the decay of the magnetization takes place ap- 
pears to be already visible within the imaginary- 
voltage data shown in Fig. [5b] Apparently, this 
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is due to the rather weak voltage-dependence of 
imaginary- voltage data (cf. figure |4aJ) . Voltages 
above lOfc^TV were not accessible to the Max- 
Ent, due to a strong sign problem occurring for 
the QMC simulations of the effective-equilibrium 
systems associated to the high-<y9 m tails. 

We again may compare the voltage dependence 
of the stationary non-equilibrium magnetization 
to the temperature dependence in equilibrium. 
Since we here are at a finite temperature T = 
Tk/2, hence the magnetization is smaller than the 
value at T = 0, the natural thing to look at is the 
curve M(T) • [M(T K /2)/M(0)} and rescale tem- 
perature with an appropriate factor. The result 
is shown as dot-dashed line in Fig. [5] Although 
one can reach a reasonable match for low voltages, 
a significant deviation occurs already at moderate 
bias. Thus there does not seem to exist a simple 
mapping $ — > T which will bring the curves to 
overlap, .i.e. it again seems doubtful that one can 
describe the effect of finite bias voltage by an ef- 
fective temperature scale, at least beyond small 
bias voltages of the order of the Kondo scale. 

On the other hand, a rather good account for 
all data can be achieved by the very simple ansatz 

m($) a 1 
B ~ B I- 2 | ~ 

V / h 2 + 4 2 

where <l := $/(2T K ). The result of this fit with 
a = 0.52, !i«2 and c ps 3 is shown as the double- 
dot-dash curve in Fig. [5] Note that this formula 
gives the right behavior in the two limits $ — >• 0, 
viz M/B oc 1 - c$ 2 with some numerical constant 
c, and $ — > oo, viz M/B oc 1/$. From scal- 
ing analysis^ one would expect that, in particular 
for large bias, additional logarithmic corrections 
appear. Due to the limited data space available 
we are of course not able to resolve those; fur- 
thermore, it is not clear if these logarithmic cor- 
rections will actually be visible in the intermedi- 
ate coupling regime studied here, due to residual 
charge fluctuations. We therefore view the above 
formula as a reasonable description in the regime 
of bias, temperature and field of the order of the 
Kondo temperature for the intermediate coupling 
regime of the SI AM. 

V. SUMMARY 

The present paper presents a detailed study on 
how the imaginary- voltage formalism proposed in 
Ref. [32] relates to Keldysh theory. Using series 
resummations, we are able to show up to all or- 
ders that static expectation values of observables, 



which satisfy certain symmetry relations with re- 
spect to the Keldysh contour, map exactly onto 
the corresponding expressions in Keldysh pertur- 
bation theory. In particular, it was pointed out 
that in order to obtain a physical expectation 
value, the limiting process iip m — > $ has to be 
taken as principal-value. This prescription en- 
sures, that one generates the principal- value in- 
tegrals which emerge in the proper real-time the- 
ory. For dynamical correlation functions, this was 
shown explicitly up to fourth order of perturba- 
tion theory 

As one important novel result of the present 
paper we were able to provide an exact spec- 
tral representation for static expectation values 
similar to a Lchmann representation. Based on 
the representation, using unbiased numerical data 
from continuous-time quantum Monte-Carlo sim- 
ulations, we found that the evaluation of the lim- 
iting procedure as principal-value expression does 
indeed give real numbers as physical expectation 
values. Consequently, the theory is found to be 
fully consistent in this respect beyond the pertur- 
bation arguments given. The double occupancy as 
function of bias voltage computed this way shows 
features similar to straight-forward second-order 
perturbation theory, but we find them to be more 
pronounced. For the magnetic susceptibility we 
were able to give numerical estimates on the de- 
struction of the Kondo effect. A comparison to 
equilibrium NRG shows that the dependence on 
bias voltage for both, the double occupancy and 
the magnetic susceptibility, cannot be explained 
by a simple effective-temperature interpretation. 
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Appendix A: Cancellation of overlapping 
(5-functions in Eq. (1511) 

With a set of {^fc^V^fcia^ = 1, ••■ ,6} 
appearing for the matrix elements in Eq. (|51l) . 
we categorize the thermal factor e^ $Y °{"' m ! > as 
follows, (i) If Y 0n = Y 0m = Yqi, Eq. (|5TJ) van- 
ishes, (ii) If only one of Y 0n7 Y 0m7 Yoi is differ- 
ent from others, (Y n,Y Qm ,Y i) g {(Y ,Y ,Y + 

1) ,(F + 1,Y ,Y ),{Y ,Y + 1,Y ),{Y ,Y ,Y + 

2) , {Y + 2, Y , Y ), (Y ,Y + 2, Y )} for some ref- 
erence value Yq. If we take the case of 
(Y 0n ,Y 0m ,Yoi) = (Yo,Y ,Yq + 1), the terms 
contributing for the matrix elements V nm , 
V ml and A ln are from ^^V^V^, 

^R kl ^Lk 2 i>l 5h ip ae ~ k6 , and 4 fca Vi? fcl ^ ^ , 

respectively, where (&!,••• , fcg) is a some per- 
mutation of (^3, /C3, ^4, ^4, • • • , fc6j &6)- The reser- 
voir indices should be chosen such that &s = a§ 
and a? = fig, and &2, &3, 014) should satisfy 
Yon = Yo m . The fii indices are summed over for 
L/R. Then the term in Eq. (|51[) becomes propor- 
tional to 

(t L t R m+tir n ifl(*i)iv w »(i-2+e^). 

i=l,6 

I 



For other combinations of (Ion, ^omj ^bz) = 
(y , % + l,Yo), (Y + 1,Y ,Y ) the thermal fac- 
tor becomes (1 - 2e' 3 * + 1) and (e' 9 * - 2 + 

1) , respectively, and all three contributions sum 
up to zero. With the case of (Y ,Y n ,Y + 

2) , the contribution becomes (Wi?) 4 ^! + 
4) 2 n i= i,6l5(^)| 2 e^ y °(l-2 + e 2 ^). The other 
terms have factors of ( 1 - 2e 2/3<1, + 1 ) , (e 2/3 * - 2 + 1 ) , 
and these sum up to zero again. 

(iii) When all of Yo n ,Yo m ,Yoi are different, 
(T n , Y 0m , Y i ) is a permutation of (Y , Y + 1 , Y + 
2). Since V, A arc at most two-particle oper- 
ators the difference of F-values between states 
cannot be greater than two. If (Yo n , Yo m , Yqi) = 
(Y ,Y + l,Y + 2), the factor in Eq. (JSTJ becomes 
proportional to 



{tLt R )\tl+t\f J] \9(h)\ 2 e^(l-2e^+e^). 

t=l,6 



Permuting (Yq, Yq + 1, Yo + 2) the sum of the ther- 
mal factors can be easily shown to be zero. 



Appendix B: Fourth order expansion of electron self-energy 



We investigate the energy-pole structure in the real-time perturbation expansion to verify that the 
^-function residue disappears and the energy denominators can be interpreted as principal-valued. 
In the following we consider the perturbation expansion for the self-energy in the fourth order of 
Coulomb parameter U, E^ (t, 0) according to the time-orderings along the Keldysh contour, FIG. [^a- 
d). Different types of time-orderings will be considered shortly. These time-orderings have one of the 
intermediate time (marked as cross) within a finite time-interval fixed by time at and t. Given a time- 
ordering, a particular Wick's contraction should be chosen. The chosen Wick's contraction is according 
to the diagrams in (g-h) which correspond to the most non-trivial vertex correction. 

We can evaluate each contribution as follows. 
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-<H — £5 +«6-i'?)si— *(-«2 +S3 +£4 — 67)^2 — i(cb 


-C6+C7)' 


(B3) 


s d 


= /1/2/3/4/5/6/7 


p — 00 

1 ds2 


/ ds\e 
Jo 


-i(«2 


— £3 — e4+«7)si-i( — ei+C4+£5— ee— in)s2-i(ei 


—62+63)* 


(B4) 



In these shorthand notation (as discussed in the main text), we omitted the expression 
u4 [I\i J d £i\9d(ei)\ 2 } which is common to all S t terms. = [1 + e^-""-*/ 2 )]- 1 and fi = l- f l . 
After some algebra, we get 

S a + S d = 2/1/2/3/4/5/6/7 j e ^i(_ e2+e3+e4+e5 _e 6 )t _ e -i( e5 -e 6 +e 7 )tj ^gg-j 

(-e 2 + e 3 + e 4 - e 7 )(ei - e 4 - e 5 + £6) 
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(g) ft) 

FIG. 9: (a-d) Real-time Keldysh contour for self-energy E^(t,0) in the fourth order perturbation when one 
intermediate time is in the finite interval [0, t] and the other time in along the contour stretching to — oo. The 
Wick's contraction is taken as shown in (g-h). The dummy label of (g) is used for time-orderings (a,c,e) and 
(h) used for (b,d,f). The cross represents the intermediate times Si and S2 for interaction, in addition to the 

creation/annihilation points and t. 



The exponential terms cancel each other at the energy poles and (£2 — £3~ £4+67) 1 and (ei— £4 — 
give well-defined principal-valued integral. This is a typical behavior since an integral within a finite 
interval (0, t) does not need the convergence factor e nt and, accordingly, principal- valued integral is 
enough. The same can be said for the combination Sb + S c . 

Now, we take the imaginary-time contours in FIG.|H]Je-f). After straightforward calculations, we have 
(ei = e. L - Oie^/2) 

p -(-e2+e3+e4+e5-«6)T _ p -(c5-e6+«7)"r 

s e = f^hhhhhjz — = — - _ L - v — (B6) 

(ei - £4 - £5 + £6j(-£2 + £3 + £4 - £7) 

p -(ei-e2+«3)T _ p — — ?4+?7)t 

-hhhhhhh T . — : — ~ ^~ u ~ ^~ ^~ — ^T- (B7) 
(ei - e 4 - £5 + £6j(-£2 + £3 + £4 - £7) 

Here flB6j| corresponds to S a of ([BT]) and (JBTJ to S c of ([B3)l . Similarly for Sf, 

g-(ei— €4+€r)r _ g — (ei-e 2 +e3)T 

Sf = /1/2/3/4/5/6/77Z z = — , ~ w — ~ — j z — —: =-r (B8) 

(£l - £4 - £5 + £e)(-£2 + £3 + £4 - £7) 

e -(e5-ee+«7)T _ g-(-e2+e3+«4+e5-e6)r 

-/1/2/3/4/5/6/77: ; ; — , - w — z — — — — =-r- (B9) 

(£1 - £4 - £5 + £6j(-£2 + £3 + £4 - £7) 

At the energy poles at for e v irj, Sf becomes identical to S e - Similarly to the real-time diagrams, 
(—62 + £3 + £4 — £7) 1 has a well-defined principal- value integral regardless of the sign of 77. Therefore 
for diagrams S a — Sf we have correct analytic continuation of imaginary-time results to those of the 
real-time via 

' zr±-—^v( (BIO) 

£1 - £4 - £5 + £6 V £1 - £4 - £5 + £6 / 
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FIG. 10: Different time-ordering with two intermediate interaction events extend to infinity. (a,d,e,g) use the 

iabei in FIG.l(g) and (b,c,f,h) FIG.l(h). 



In FIG. 1101 we consider the remaining time-orderings with the two intermediate interaction points 
extending to infinity. These are harder to deal with, as we discuss below, since the energy poles may 
overlap. 



D a = hhhhhhh // dsids2e-^ ei - e *- e5 + £6 -^ sl - i ^- £2+e3 - e5 + £6 - £7 -^ S2 - i ( £5 - e6+£7 )*(Bii) 



D, 
D d 
A 
Df 



— /1/2/3/4/5/6/7 / e _ ^ _Cl+e2_e3+e5_e6+e7_il '^ 1_ ^ _ei+e4+e5_e6_i ''^ S2_ ^ Cl_e2+C3 ^ (B12) 

fi /2/3 Ja fhf&h J e _i ^ ei+e2_e3+e5_e6+£7_i?, ^ 1_i ^ ei+e4+C5_C6_i?, ^ 2_ ^ e5 ~ e6+e7 - ) * (B13) 

— /1/2/3/4/5/6/7 J e _ ^ ei_e4_e5+e6_ ^^ 1_ ^ ei_e2+e3_e5+e6_C7_!;,? ' S2_ ^ ei_e2+e3 ^ (B14) 

—fxfihfifbhh J e~^ ei_C4_e5+e6_i?, ^ 1_ ^ _e2+e3+e4_e7_ "'' S2_ ^ ei_e4+e7 ^ (B15) 

fi /2/3 Jifsfsh / e _ ^ C2_e3_e4+£7_i, '^ 1_ ^ _ei+C4+e5_e6_iT '^ 2_i ' _e2+e3+C4+£5_e6 '* (B16) 



After integrals over s\ and S2 it is easy to see that D a (irj) = D c (—irj) and D^irj) = D^—irj). For D e 
and Df, we can swap the dummy indices as 1 -s-> 7, 2 <-> 6, and 3 <-> 5, and it becomes D e (iri) = D e (—irj) 
and Df(irj) = Df(—irj). Therefore, we obtain the desired result as (|B10[) . 



J2 D k (ir 1 )=J2D k (-ir 1 )=J2VD k (±iv)- 

k—a,- ■• ,f k k 



(B17) 



In deriving these relations, no assumptions of L/R and particle- hole symmetry have been used. One 
can rewrite D a as 



D a — /1/2/3/4/5/6/7 



e -«(e5-e6+e7)t 
£2 — £3 — e 4 + £7 



ei - e 2 + e 3 - e 5 + e 6 - e 7 - 177 ei - e 4 - e 5 + e 6 - z?/ 



. (B18) 

Here the +irj in the denominator will be cancelled by D c and all fractions can be written as principal- 
valued, unless the poles coincide. 
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We can now turn to the imaginary-time diagrams FIG. [TUJg,h). 

D _ hhhfihhh ( 1 . 1 \ e _(g 5 _ g6+ e 7 ) T 

3 -(?2 - h - e 4 + e» \ h - h + h - h + h - h ei - e 4 - e 5 + e 6 / 

hfvhhhhh e -(gi-.- a +.-3)r ^ hhhhhhh e-^-^+j^ 

(e 2 -h~h + h) (h - h + h ~ h + h - h) {h - e 3 - e 4 + e 7 ) (?i - e 4 - e 5 + e 6 j 

After swapping 1 f> 7. 2 f> 6. and 3 O 5, the first two terms correspond to D a and Z? c for — s- 
177 an< i the third term to D e . Using a similar technique in (|B18|) . we can decouple the product of 
energy denominators to a sum of simple poles of e v and then by taking the limit Eq. (|B10[) , all energy 
denominators become principal- valued, unless poles coincide. 

Now we deal with the case when the 8- functions overlap. As discussed in section UlI A[ the double-<5 
terms manifest as terms proportional to T 2 . The terms D a , D c and D e have double-5 terms cancelled 
among themselves. At the energy-poles e± — e 4 — £5 + eg = and £2 — £3 — e 4 + e- = 0, 

D a = d c oc fihhhhhf7Y e ~ i{e5 ~ e6+eT)t - ( B2 °) 



For D e , we first rewrite 



dsi = / dsi + ds u (B21) 



and note that the second integral with a finite interval should not contribute a 5-function. So as long 
as double-<5 is concerned, we only consider the first interval, 

D e cx -fihhUhf 6 f7T 2 e- i ^- e ° + ^ t -> -hhhhhhhT 2 e- l[ ^ B+e7)t . (B22) 

where at the last step the dummy indices are swapped as 1 -f-> 5 and 4 f> 6. Therefore, the double-5 
terms disappear in D a + D c + D e . The same is true with + Dd + Df, and it shows that the all energy 
poles for the fourth-order vertex corrections, FIG. [9]^ g-h), are interpreted as principal- valued. 
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